On cherche à étudier l’effet de trois facteurs sur le transcriptome des racines d’Arabidopsis thaliana et de la micro Tomate.
load(paste0("GenesCO2_", specie, ".RData"))
# quantification file
data <- read.csv("quantifFiles/QuantifGenes.csv", h = T, sep = ",")
rownames(data) <- data$Gene
genes = which(!(grepl("__", rownames(data))))
not_quant = data[which((grepl("__", rownames(data)))), ]
data = data[genes, grepl("R", colnames(data))]
keep <- rowSums(data) >= 10
data <- data[keep, ]
group <- sapply(colnames(data), getLabel, with.rep = F)
colnames(data) <- sapply(colnames(data), getLabel)
specie = "At"
clusteredGenes <- clustering(sharedBy3, data)****************************************
coseq analysis: Poisson approach & none transformation
K = 2 to 12
Use set.seed() prior to running coseq for reproducible results.
****************************************
Running g = 2 ...
[1] "Initialization: 1"
[1] "Log-like diff: 0"
Running g = 3 ...
[1] "Initialization: 1"
[1] "Log-like diff: 0"
Running g = 4 ...
[1] "Initialization: 1"
[1] "Log-like diff: 9.85664883046411e-11"
Running g = 5 ...
[1] "Initialization: 1"
[1] "Log-like diff: 5.6843418860808e-13"
Running g = 6 ...
[1] "Initialization: 1"
[1] "Log-like diff: 5.98125211581646e-08"
Running g = 7 ...
[1] "Initialization: 1"
[1] "Log-like diff: 0"
Running g = 8 ...
[1] "Initialization: 1"
[1] "Log-like diff: 0"
Running g = 9 ...
[1] "Initialization: 1"
[1] "Log-like diff: 3.95914412365528e-11"
Running g = 10 ...
[1] "Initialization: 1"
[1] "Log-like diff: 1.28659394249553e-09"
Running g = 11 ...
[1] "Initialization: 1"
[1] "Log-like diff: 0"
Running g = 12 ...
[1] "Initialization: 1"
[1] "Log-like diff: 0.0035332049280612"
[1] "Log-like diff: 0.000485432284904164"
[1] "Log-like diff: 6.66189752109858e-05"
[1] "Log-like diff: 9.14113201133659e-06"
$ICL
$profiles
$boxplots
$probapost_barplots
*************************************************
Model: Poisson
Transformation: none
*************************************************
Clusters fit: 2,3,4,5,6,7,8,9,10,11,12
Clusters with errors: ---
Selected number of clusters via ICL: 12
ICL of selected model: -740101.6
*************************************************
Number of clusters = 12
ICL = -740101.6
*************************************************
Cluster sizes:
Cluster 1 Cluster 2 Cluster 3 Cluster 4 Cluster 5 Cluster 6 Cluster 7
11 11 20 12 13 18 6
Cluster 8 Cluster 9 Cluster 10 Cluster 11 Cluster 12
10 9 5 14 2
Number of observations with MAP > 0.90 (% of total):
131 (100%)
Number of observations with MAP > 0.90 per cluster (% of total per cluster):
Cluster 1 Cluster 2 Cluster 3 Cluster 4 Cluster 5 Cluster 6 Cluster 7
11 11 20 12 13 18 6
(100%) (100%) (100%) (100%) (100%) (100%) (100%)
Cluster 8 Cluster 9 Cluster 10 Cluster 11 Cluster 12
10 9 5 14 2
(100%) (100%) (100%) (100%) (100%)
Class: pca dudi
Call: dudi.pca(df = log(data + 0.1), center = TRUE, scale = TRUE, scannf = FALSE,
nf = 4)
Total inertia: 24
Eigenvalues:
Ax1 Ax2 Ax3 Ax4 Ax5
17.2837 4.7082 0.6320 0.5028 0.2617
Projected inertia (%):
Ax1 Ax2 Ax3 Ax4 Ax5
72.015 19.617 2.633 2.095 1.090
Cumulative projected inertia (%):
Ax1 Ax1:2 Ax1:3 Ax1:4 Ax1:5
72.02 91.63 94.27 96.36 97.45
(Only 5 dimensions (out of 24) are shown)
NULL
load("./normalized.count_At.RData")
log.data <- log2(normalized.count[sharedBy3, ] + 1)
Norm.interest.corr <- corr.test(t(log.data), method = "pearson", ci = F)
Norm.interest.corr$p[lower.tri(Norm.interest.corr$p, diag = -TRUE)] = NA
Pval.adj <- as.data.frame(as.table(Norm.interest.corr$p))
Norm.interest.corr$r[lower.tri(Norm.interest.corr$r, diag = TRUE)] = NA
Correlation <- as.data.frame(as.table(Norm.interest.corr$r))
Cor.table <- na.exclude(cbind(Correlation, Pval.adj))[, c(1, 2, 3, 6)]
colnames(Cor.table) <- c("gene1", "gene2", "cor", "p.adj")
Cor.table.filt <- Cor.table[(abs(Cor.table[, 3]) > 0.9 & Cor.table[, 4] < 0.01),
]
g <- graph.data.frame(Cor.table.filt[, 1:2], directed = -FALSE)
V(g)$color <- clusteredGenes[V(g)]
degree <- degree(g)
hist(degree, breaks = 30)Node_nw_st <- data.frame(degree, betweenness)
plot.igraph(g, vertex.size = 5, vertex.label.cex = 0.3, color = clusteredGenes) # Nitrate
****************************************
coseq analysis: Poisson approach & none transformation
K = 2 to 12
Use set.seed() prior to running coseq for reproducible results.
****************************************
Running g = 2 ...
[1] "Initialization: 1"
[1] "Log-like diff: 4.59245086403826e-10"
Running g = 3 ...
[1] "Initialization: 1"
[1] "Log-like diff: 0.02732707907575"
[1] "Log-like diff: 0.000268834437022747"
[1] "Log-like diff: 3.19062163711692e-06"
Running g = 4 ...
[1] "Initialization: 1"
[1] "Log-like diff: 0.00317565071129877"
[1] "Log-like diff: 6.45337854905392e-05"
[1] "Log-like diff: 1.32643040373637e-06"
Running g = 5 ...
[1] "Initialization: 1"
[1] "Log-like diff: 780.379968968336"
[1] "Log-like diff: 543.197241442135"
[1] "Log-like diff: 752.743987958856"
[1] "Log-like diff: 533.723082605684"
[1] "Log-like diff: 94.16844828535"
Running g = 6 ...
[1] "Initialization: 1"
[1] "Log-like diff: 0.429716284972326"
[1] "Log-like diff: 0.0110679937471261"
[1] "Log-like diff: 0.000294696082274726"
[1] "Log-like diff: 7.96987497508894e-06"
Running g = 7 ...
[1] "Initialization: 1"
[1] "Log-like diff: 6.9358451781909e-07"
Running g = 8 ...
[1] "Initialization: 1"
[1] "Log-like diff: 1.48705296112439e-05"
[1] "Log-like diff: 4.44065548776962e-06"
Running g = 9 ...
[1] "Initialization: 1"
[1] "Log-like diff: 0.187270797612783"
[1] "Log-like diff: 0.0804456928053696"
[1] "Log-like diff: 0.0318910951551601"
[1] "Log-like diff: 0.0125871810649834"
[1] "Log-like diff: 0.00495944033215423"
Running g = 10 ...
[1] "Initialization: 1"
[1] "Log-like diff: 4.49220970679676e-05"
[1] "Log-like diff: 5.0755175600159e-06"
Running g = 11 ...
[1] "Initialization: 1"
[1] "Log-like diff: 2.50112567030669e-05"
[1] "Log-like diff: 1.44340664931519e-06"
Running g = 12 ...
[1] "Initialization: 1"
[1] "Log-like diff: 3.99837858111596e-06"
$ICL
$profiles
$boxplots
$probapost_barplots
*************************************************
Model: Poisson
Transformation: none
*************************************************
Clusters fit: 2,3,4,5,6,7,8,9,10,11,12
Clusters with errors: ---
Selected number of clusters via ICL: 12
ICL of selected model: -3013888
*************************************************
Number of clusters = 12
ICL = -3013888
*************************************************
Cluster sizes:
Cluster 1 Cluster 2 Cluster 3 Cluster 4 Cluster 5 Cluster 6 Cluster 7
109 24 105 49 122 7 121
Cluster 8 Cluster 9 Cluster 10 Cluster 11 Cluster 12
35 21 158 32 54
Number of observations with MAP > 0.90 (% of total):
836 (99.88%)
Number of observations with MAP > 0.90 per cluster (% of total per cluster):
Cluster 1 Cluster 2 Cluster 3 Cluster 4 Cluster 5 Cluster 6 Cluster 7
109 24 105 49 122 7 121
(100%) (100%) (100%) (100%) (100%) (100%) (100%)
Cluster 8 Cluster 9 Cluster 10 Cluster 11 Cluster 12
35 21 157 32 54
(100%) (100%) (99.37%) (100%) (100%)
Class: pca dudi
Call: dudi.pca(df = log(data + 0.1), center = TRUE, scale = TRUE, scannf = FALSE,
nf = 4)
Total inertia: 24
Eigenvalues:
Ax1 Ax2 Ax3 Ax4 Ax5
19.1693 3.3531 0.5281 0.3930 0.1205
Projected inertia (%):
Ax1 Ax2 Ax3 Ax4 Ax5
79.872 13.971 2.200 1.638 0.502
Cumulative projected inertia (%):
Ax1 Ax1:2 Ax1:3 Ax1:4 Ax1:5
79.87 93.84 96.04 97.68 98.18
(Only 5 dimensions (out of 24) are shown)
NULL
log.data <- log2(normalized.count[sharedBy3, ] + 1)
Norm.interest.corr <- corr.test(t(log.data), method = "pearson", ci = F)
Norm.interest.corr$p[lower.tri(Norm.interest.corr$p, diag = -TRUE)] = NA
Pval.adj <- as.data.frame(as.table(Norm.interest.corr$p))
Norm.interest.corr$r[lower.tri(Norm.interest.corr$r, diag = TRUE)] = NA
Correlation <- as.data.frame(as.table(Norm.interest.corr$r))
Cor.table <- na.exclude(cbind(Correlation, Pval.adj))[, c(1, 2, 3, 6)]
colnames(Cor.table) <- c("gene1", "gene2", "cor", "p.adj")
Cor.table.filt <- Cor.table[(abs(Cor.table[, 3]) > 0.9 & Cor.table[, 4] < 0.01),
]
g <- graph.data.frame(Cor.table.filt[, 1:2], directed = -FALSE)
V(g)$color <- clusteredGenes[V(g)]
degree <- degree(g)
hist(degree, breaks = 30)Node_nw_st <- data.frame(degree, betweenness)
plot.igraph(g, vertex.size = 5, vertex.label.cex = 0.01, color = clusteredGenes)****************************************
coseq analysis: Poisson approach & none transformation
K = 2 to 12
Use set.seed() prior to running coseq for reproducible results.
****************************************
Running g = 2 ...
[1] "Initialization: 1"
[1] "Log-like diff: 1.07117870129514e-10"
Running g = 3 ...
[1] "Initialization: 1"
[1] "Log-like diff: 6.38532487400312e-05"
[1] "Log-like diff: 6.09656778216561e-06"
Running g = 4 ...
[1] "Initialization: 1"
[1] "Log-like diff: 0.0111171821348677"
[1] "Log-like diff: 0.00141141140834478"
[1] "Log-like diff: 0.000161399031195941"
[1] "Log-like diff: 2.23980092322051e-05"
[1] "Log-like diff: 2.78977087120325e-06"
Running g = 5 ...
[1] "Initialization: 1"
[1] "Log-like diff: 3567.60751860177"
[1] "Log-like diff: 3726.79928110929"
[1] "Log-like diff: 1847.6149884063"
[1] "Log-like diff: 8431.18437529454"
[1] "Log-like diff: 2533.24931846412"
Running g = 6 ...
[1] "Initialization: 1"
[1] "Log-like diff: 61.6607030025398"
[1] "Log-like diff: 220.279999620299"
[1] "Log-like diff: 343.597475596329"
[1] "Log-like diff: 118.657798253033"
[1] "Log-like diff: 89.3377189676472"
Running g = 7 ...
[1] "Initialization: 1"
[1] "Log-like diff: 294.715502575396"
[1] "Log-like diff: 150.39616865371"
[1] "Log-like diff: 197.001670751701"
[1] "Log-like diff: 20.6833498346011"
[1] "Log-like diff: 6.69058247831302"
Running g = 8 ...
[1] "Initialization: 1"
[1] "Log-like diff: 43.5314394208409"
[1] "Log-like diff: 327.586113284712"
[1] "Log-like diff: 500.912263648639"
[1] "Log-like diff: 156.206329141795"
[1] "Log-like diff: 184.589709254823"
Running g = 9 ...
[1] "Initialization: 1"
[1] "Log-like diff: 481.077552880626"
[1] "Log-like diff: 58.6942035806858"
[1] "Log-like diff: 992.988620924576"
[1] "Log-like diff: 30.1221236683701"
[1] "Log-like diff: 235.436537028793"
Running g = 10 ...
[1] "Initialization: 1"
[1] "Log-like diff: 477.227841689411"
[1] "Log-like diff: 126.914316995362"
[1] "Log-like diff: 516.989617982976"
[1] "Log-like diff: 211.394809733677"
[1] "Log-like diff: 38.5426970431339"
Running g = 11 ...
[1] "Initialization: 1"
[1] "Log-like diff: 1.94951214516996"
[1] "Log-like diff: 185.070654937831"
[1] "Log-like diff: 28.381146970325"
[1] "Log-like diff: 26.0836450267923"
[1] "Log-like diff: 44.8672622161952"
Running g = 12 ...
[1] "Initialization: 1"
[1] "Log-like diff: 145.645179740505"
[1] "Log-like diff: 37.4337944851441"
[1] "Log-like diff: 391.669701764039"
[1] "Log-like diff: 396.307727963539"
[1] "Log-like diff: 300.789765029271"
$ICL
$profiles
$boxplots
$probapost_barplots
*************************************************
Model: Poisson
Transformation: none
*************************************************
Clusters fit: 2,3,4,5,6,7,8,9,10,11,12
Clusters with errors: ---
Selected number of clusters via ICL: 12
ICL of selected model: -3391000
*************************************************
Number of clusters = 12
ICL = -3391000
*************************************************
Cluster sizes:
Cluster 1 Cluster 2 Cluster 3 Cluster 4 Cluster 5 Cluster 6 Cluster 7
688 645 110 220 133 30 73
Cluster 8 Cluster 9 Cluster 10 Cluster 11 Cluster 12
303 122 110 366 41
Number of observations with MAP > 0.90 (% of total):
2773 (97.61%)
Number of observations with MAP > 0.90 per cluster (% of total per cluster):
Cluster 1 Cluster 2 Cluster 3 Cluster 4 Cluster 5 Cluster 6 Cluster 7
682 638 105 215 122 28 73
(99.13%) (98.91%) (95.45%) (97.73%) (91.73%) (93.33%) (100%)
Cluster 8 Cluster 9 Cluster 10 Cluster 11 Cluster 12
291 118 107 354 40
(96.04%) (96.72%) (97.27%) (96.72%) (97.56%)
Class: pca dudi
Call: dudi.pca(df = log(data + 0.1), center = TRUE, scale = TRUE, scannf = FALSE,
nf = 4)
Total inertia: 24
Eigenvalues:
Ax1 Ax2 Ax3 Ax4 Ax5
22.03089 1.11520 0.27107 0.11321 0.06945
Projected inertia (%):
Ax1 Ax2 Ax3 Ax4 Ax5
91.7954 4.6467 1.1295 0.4717 0.2894
Cumulative projected inertia (%):
Ax1 Ax1:2 Ax1:3 Ax1:4 Ax1:5
91.80 96.44 97.57 98.04 98.33
(Only 5 dimensions (out of 24) are shown)
NULL
log.data <- log2(normalized.count[sharedBy3, ] + 1)
Norm.interest.corr <- corr.test(t(log.data), method = "pearson", ci = F)
Norm.interest.corr$p[lower.tri(Norm.interest.corr$p, diag = -TRUE)] = NA
Pval.adj <- as.data.frame(as.table(Norm.interest.corr$p))
Norm.interest.corr$r[lower.tri(Norm.interest.corr$r, diag = TRUE)] = NA
Correlation <- as.data.frame(as.table(Norm.interest.corr$r))
Cor.table <- na.exclude(cbind(Correlation, Pval.adj))[, c(1, 2, 3, 6)]
colnames(Cor.table) <- c("gene1", "gene2", "cor", "p.adj")
Cor.table.filt <- Cor.table[(abs(Cor.table[, 3]) > 0.9 & Cor.table[, 4] < 0.01),
]
g <- graph.data.frame(Cor.table.filt[, 1:2], directed = -FALSE)
V(g)$color <- clusteredGenes[V(g)]
degree <- degree(g)
hist(degree, breaks = 30)Node_nw_st <- data.frame(degree, betweenness)
plot.igraph(g, vertex.size = 5, vertex.label.cex = 0.01, color = clusteredGenes)